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ABSTRACT 


A  dynamical-numerical  experiment  is  designed  to  evolve  thermally- 
driven  circulations,  symmetric  about  the  polar  axis,  in  the  region  from 
equator  to  pole.  The  selection  of  variables  and  transformations  of  the 
hydrostatic  equation  system  (primitive  equations)  and  the  development 
of  the  numerical  analogue  are  discussed  in  some  detail.  The  region  is 
discretized  into  volume  modules  by  nineteen  levels  in  the  vertical  and 
by  forty-nine  equal  intervals  from  equator  to  pole.  The  upper  boundary 
is  essentially  unrestricted.  Mass  and  absolute  angular  momentum  are 
treated  conservatively  and  in  terms  of  fluxes  across  module  interfaces. 
Heating  rate  distributions  are  prescribed  in  space  and  are  stepwise 
introduced  to  disturb  an  atmosphere  at  relative  rest  or  in  balanced 
zonal  motion.  The  fields  are  portrayed  in  a  representative  meridional 
plane.  Pertinent  numerical  analyses  are  included  and  compressior- 
gravity  oscillations  are  analyzed. 

The  first  experiment,  designed  as  a  test,  evolves  the  circulation 
produced  by  uniformly  heating  a  narrow  latitude  zone  at  mid-latitude  in 
a  Standard  Atmosphere  at  relative  rest.  Although  receiving  thermal 
energy,  the  heated  zone  locally  cools  due  to  the  combined  effect  of  the 
resulting  circulation. 

In  the  major  experiment  a  heating  and  cooling  distribution,  approxi¬ 
mating  mean  winter  conditions,  disturbs  a  Standard  Atmosphere  at  relative 
rest.  The  integration  evolves  the  initiation  of  a  planetary  circulation 
as  frequently  discussed  qualitatively  in  elementary,  text  books.  Various 
applications  of  the  model  and  further  improvements  are  indicated. 
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I  INTRODUCTION 


The  addition  and  renoval  of  thermal  energy  in  the  atmosphere  (i.e., 
non -adiabatic  heating)  is  exceedingly  complex,  depending  as  it  does  on 
radiational  exchange  processes,  convection  from  the  surface,  and  the 
storage  and  release  of  latent  heat,  and  on  other  processes  which  are 
significant  in  the  high  atmosphere.  The  heating  itself  is  for  the  most 
part  strongly  dependent  on  the  evolving  state  of  the  atmosphere,  on  the 
distribution  of  radiating  and  absorbing  substances,  the  temperature,  the 
convective  stability,  the  water  forms,  the  vertical  motion,  etc.,  and 
on  the  time  of  day.  However,  before  attempting  to  introduce  any  of  this 
complexity  into  dynamical-numerical  models,  it  is  necessary  to  have  some 
appreciation  of  both  (1)  the  distributions  of  heating  and  cooling  in 
the  atmosphere  in  terms  of  magnitude  and  location  relative  to  circula¬ 
tion  systems  and  geography,  and  (2)  the  dynamical  significance  of  such 
distributions.  The  present  project  set  out  to  contribute  to  knowledge 

y  g)  ^ 

in  both  of  these  areas.  Scientific  Reports  1  and  2  on  this  contract 
deal  with  the  former;  and  this  work  is  continuing.  The  present  report 
(Scientific  Report  3)  pertains  to  the  latter. 

The  experiments  discussed  herein  were  begun  for  the  purpose  of 
investigating  the  use  of  a  suitable  dynamical-numerical  model  for  the 
incorporation  of  heating,  and  to  learn  about  the  dynamical  significance 
of  heating  distributions.  In  these  experiments  the  heating  field  is 
independently  prescribed,  by  choice,  in  space  and  time  that  is,  the 
heating  will  not  be  coupled  to  any  of  the  evolving  parameters. 

It  was  desired  that  these  investigations  should  relate  and  con¬ 
tribute,  in  so  far  as  possible,  to  techniques  in  the  general  compre¬ 
hensive  dynamical-numerical  modeling  of  atmospheric  evolutions. 


*  References  are  listed  at  the  end  of  the  report. 
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In  this  work  the  attitude  is  to  approxinate  the  fields  as  continua.  Re¬ 
strictions  which  limit  the  continuum  are  imposed  only  in  order  to  limit 
the  magnitude  of  the  task  while  still  maintaining  enough  degrees  of 
freedom  for  certain  investigations. 

The  pertinent  laws  of  physics,  modified  by  the  hydrostatic  relation¬ 
ship,  lead  to  the  so-called  primitive  (or  hydrostatic)  equation  system. 
The  system  is  relatively  unrestricted  and  provides  a  broad  physical 
framework  for  the  progressive  incorporation  of  additional  pertinent 
processes  and  influences.  In  the  last  few  years  it  has  been  numerically 
integrated  and  investigated  by  several  groups  and  researchers.  The 
choice  of  dependent  and  independent  variables  leads  to  various  transfor¬ 
mations  of  the  system.  Our  choice  is  to  use  density  as  a  function  of 
geometric  height  as  the  principal  dependent  variable  of  state.  The 
preference  for  this  choice  and  the  system  were  discussed  in  an  earlier 
(1960)  work*.3*  In  the  system  the  heating  rate  appears  in  Richardson’s 
equation,  which  diagnostic  equation  specifies  the  field  of  vertical 
motion. 

The  dynamical  significance  of  heating  in  the  atmosphere  ranges 
through  all  scales.  Understanding  is  needed  of  its  effectiveness  in 
(1)  creating  thermally  driven  circulations,  (2)  modifying  the  evolution 
of  circulation  systems,  and  (3)  triggering  dynamic  instabilities. 

One  of  the  major  problems  in  the  use  of  the  hydrostatic  system  is 
the  determination  of  initial  numerical  data  for  the  principal  dependent 
variables  of  density,  and  horizontal  velocity.  A  smooth  specification 
of  these  fields,  in  the  scale  of  the  discretization,  is,  in  general, 
not  sufficient  because  too  much  energy  appears  in  compression-gravity 
waves,  particularly  in  the  shorter  components.  The  compression-gravity 
mechanism  plays  a  major  role  in  the  evolution.  The  problem  focuses  on 
having  the  correct  field  of  dynamic  imbalance,  as  calculated  by  the 
numerical  analogue  of  the  equation  system. 

In  designing  experiments  for  specific  investigations  the  objective 
is  a  minimum  amount  of  complication  for  a  maximum  return  of  knowledge 
within  a  specified  task  magnitude.  Along  this  theme,  certain 
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coaplications  were  selected  and  incorporated  by  design  while  other 
coaplications  were  tentatively  set  aside.  Experiments  can  always  be 
iaproved  upon,  especially  afterwards.  Later  sections  of  this  report 
contain  criticises  and  recoaaendations. 

It  was  decided  to  forego  variability  in  one  horizontal  dimension  in 
favor  of  aany  levels  in  the  vertical.  The  meridional  plane  is  taken  as 
the  field  of  variability,  and  symmetry  about  the  polar  axis  is  imposed. 
The  model  extends  from  pole  to  equator  with  no  flow  across  the  equatorial 
wall.  This  calls  for  the  use  of  spherical  coordinates. 

The  problem  of  initial-data  specification  was  simplified  by  limiting 
ourselves  to  cases  which  are  initially  in  static  (rest)  or  dynamic 
(steady)  equilibrium.  The  equilibrium  is  then  upset  by  a  specified 
distribution  of  heating.  The  numerical  analogue  permits  ready  determi¬ 
nation  of  a  field  of  zonal  motion  which  is  in  balance  with  an  arbitrarily 
specified  density  distribution.  The  integration  itself  then  leads  to 
more  general  states  for  the  study  of  subsequent  evolutions  under  the 
influence  of  heating. 

Frictional  processes,  both  internal  and  boundary,  were  also  set 
aside.  Friction  is  not  an  inseparable  complement  to  heating,  but  is  an 
additional  process  acting  in  the  evolution,  and  may  be  independently 
considered.  It  may  be  argued  that  without  boundary  friction,  there  is 
no  cognizance  in  the  model  of  the  rotation  of  the  underlying  smooth 
earth.  However,  as  designed,  the  model  conserves  total  angular  momentum 
about  the  polar  axis  at  the  value  initially  specified.  In  this  way  the 
rotation  of  the  system  is  built  into  the  evolution.  The  model  is,  of 
course,  limited  without  the  incorporation  of  friction. 

The  numerical  analogue  was  designed  in  terms  of  volume  modules, 

budgets,  and  fluxes  for  mass  and  angular  momentum.  In  this  way  total 

mass  and  total  angular  momentum  are  conserved.  Truncation  errors  in 

mass  and  angular  momentum  tendencies  cancel  internally  in  any  sub-region. 

(4) 

This  is  particularly  useful  in  pressure-tendency  calculations. 

The  conventional  "second-order-centered"  method  of  time  integration 
(i.e.,  the  leap-frog  method)  was  used  with  reservations,  but  as  it 
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turned  out,  with  aoae  advantage  in  the  study  of  the  shortest-period 
c depression-gravity  node. 

The  nodel  is  designed  to  evolve,  with  limitations,  axially  symmetric 
planetary  circulations  under  the  influence  of  specified  n on-adiabatic 
heating  distributions  which  vary  with  latitude,  altitude,  and  in  tine. 

For  Investigations  of  soae  effects  of  heating  at  high  levels,  such  as 
conceived  in  solar-weather  relationship  studies,  the  upper  boundary 
conditions  are  of  special  interest.  They  are  essentially  non-restrictive 
and  complete.  Thermal  energy  may  be  added  to,  or  removed  from,  the 
upper  reaches  of  the  atmosphere  in  detail  as  limited  by  the 
discretization. ^ 

For  the  numerical  analogue,  the  meridional  plane  is  discretized 
into  volume  modules  by  49  equal  intervals  from  equator  to  pole  and  by 
19  selected  levels  with  irregular  spacing  in  the  vertical.  A  time  step 
of  3  minutes  is  used  in  the  final  integrations.  The  heating  is  speci¬ 
fied  for  each  of  the  volume  modules. 

In  the  Integrations,  abrupt  changes  in  heating  are  avoided  by 
introducing  the  changes  gradually  over  several  time  steps.  Only  1/32  of 
the  change  is  introduced  at  first.  It  is  then  progressively  doubled 
with  successive  time  steps  until  the  new  rate  is  achieved.  This  pro¬ 
cedure  serves  two  purposes.  It  reduces  the  compression-gravity  dis¬ 
turbances  which  would  arise  by  the  larger  abrupt  changes.  It  is  also 
required,  especially  in  beginning  the  integration  from  equilibrium, 
in  order  to  minimize  arousing  the  extraneous  degrees  of  freedom  which 
enter  by  use  of  the  "second-order-centered"  integrating  procedure/6^ 

In  proceeding  from  time  zero  to  a  steady  heating  rate,  the  full  heating 
may  be  considered  to  have  begun,  effectively,  after  an  interval  of 
about  4-1/2  time-steps. 

Only  a  few  experiments  were  carried  out  with  the  model  as  developed 
and  programmed  for  a  Burroughs  220  electronic  computer.  More  was  not 
done  because  it  was  felt  that  an  optimum  had  been  reached  in  considera¬ 
tion  of  experimental  findings,  the  desirability  of  model  extensions, 
and  the  economic  advantages  of  reprogramming  for  a  larger  computer. 
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The  first  set  of  experiments  was  designed  primarily  to  test  and 
study  the  model.  Initially  the  atmosphere  is  at  rest  relative  to  the 
rotating  earth.  One  column  (i . e. ,  the  region  between  two  latitude  walls , 
not  quite  2°  apart)  at  aid- latitude  is  heated  at  a  steady  rate  which, 
were  the  column  confined  by  latitude  walls  so  that  each  level  expanded 
at  constant  pressure,  would  warm  the  column  by  10 °C  in  24  hours.  The 
general  results  bear  on  computational  stability  and  the  nature  of  the 
shortest-period  compression-gravity  mode.  The  particular  results  of 
the  experiment  showed  the  somewhat  surprising  result  that,  as  a  conse¬ 
quence  of  the  thermally  forced  circulation,  the  heated  column  becomes 
colder  than  the  unheated  environment.  This  result  is  not  in  error  and 
is  physically  consistent.  It  demonstrates  dramatically  the  dangers  of 
jumping  to  conclusions  by  qualitative  reasoning. 

The  major  experiment  also  was  begun  with  an  atmosphere  at  rest 
relative  to  the  rotating  earth.  The  motion  is  then  forced  by  a  steady 
heating  and  cooling  distribution  patterned  after  the  winter  season. 

The  results  are  suitable  as  a  demonstration  of  common  elementary  text¬ 
book  discussions  of  the  atmosphere  as  a  heat  engine  and  the  formation 
of  zones  of  easterly  and  westerly  winds. 

The  circulations  developed  by  the  model  must  be  interpreted  with 
the  realization  that  these  circulations  continue  to  develop  without  being 
permitted  to  break  down  into  longitudinal  disturbances  due  to  dynamic 
instability.  The  zonal  circulation  and  baroclinicity  fields  progress 
beyond  thresholds  of  such  instability. 

If  the  temptation  arises  to  compare  results  with  mean  winter-season 
meridional  fields,  it  should  be  kept  in  mind  that  the  real  data  contain 
the  effects  of  longitudinal  waves  and  eddies.  Furthermore,  the  model 
is  a  limited  dynamical-numerical  approximation.  Lastly,  the  integration 
has  been  carried  out  only  far  enough  to  show  a  few  derivatives  in  the 
sense  of  a  Taylor  Series  time  expansion. 

The  model  is  general  in  the  sense  that  it  can  be  applied  to 
arbitrary  planets  by  changing  dimensions,  parameters,  and  initial  values. 
Such  integrations  would  give,  under  specified  thermal  forcing  fields, 
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soae  preliminary  notions  of  the  planet's  circulation.  The  regions  of 
developing  shear  and  baroelinicity  may  give  soae  indication  of  cyclone 
belts. 

Another  application  is  that  the  model  can  be  used  for  saaller-scale 
axially  symmetric  circulations  Including  the  convective  scale.  This  is 
done  by  shrinking  the  horizontal  dimension  and  imposing  a  constant 
Coriolis  force.  In  this  way  the  effects  of  a  strong  local  release  of 
thermal  energy  might  be  examined.  Application  is  limited  by  the 
imposition  of  axial  symmetry. 
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II  THE  MATHEMATICAL-PHYSICAL  SYSTEM 


The  basic  physical  relationships  are  expressed  by  the  equation  tor 
continuity  of  natter, 

v  •  (PV)  (1) 

the  acceleration  equation, 

.  y.  v  v  -  2  «x  V  -  ivp-  v*  (2) 

and  the  thermodynamic  equation, 

j=K  In  p p~K  =  ^  $  (3) 

The  symbols  have  their  usual  meaning:  p  is  the  density,  y  is  the  three- 
dimensional  velocity  vector,  8  is  the  earth's  rotation  vector,  p  is 
the  pressure,  ♦  is  the  apparent  gravity  potential,  Cy  is  the  specific 
heat  of  air  at  constant  volume,  K  -  Cp/Cv  where  is  the  specific 

heat  of  air  at  constant  pressure,  R  is  the  gas  constant  for  dry  air, 
and  Q  is  the  rate  at  which  thermal  energy  is  being  added.  Equation  (3) 
also  implicitly  contains  the  perfect  gas  law, 

p  =  p  RT 

where  T  is  the  absolute  temperature. 

As  spherical  coordinates,  we  take  <p  as  the  latitude,  X  the  longi¬ 
tude,  and  r  the  radial  distance  from  the  center  point.  The  corresponding 


unit  vectors  are 


i  =  r  cos  <f>  VX 
j  =  r  V  <*> 

V  =  Vr 


The  velocity  vector  and  the  del  operator,  in  corresponding  components, 
are 


V 

V 


ui  +  v  $  +  wk 

_ i _  i  *ii  +k  »-  . 

r  cos  <p  dX  r  8 <f>  dr 


The  expansion  of  the  system  into  spherical  coordinates  is  readily  ac¬ 
complished  by  first  expanding  the  vector  pair  (or  dyadic)  VV  ,  which 
we  shall  exhibit  here,  in  matrix  form,  for  reference: 


i 

1  8u  v  . 

<p  +  -  | 
^  r  1 

i 

1 

i 

9v  A 

-  tan  <p 
r 

Ic 

11  9w  u  ’ 

r  cos  <(>  8A  r 

r  cos  <p  8X 

i  r  cos  <f>  dX  r 

i 

1  9u 

i 

1 

8v 

w 

i 

i  I  _  I 

r  9  <p 

i 

i 

r 

8  <f> 

+  “ 
r 

r  dip  r 

i 

qjIqj 
►1  Ic 

i 

1 

1 

8v 

8r 

i 

1  dw 

i  dr 

The  earth  is  taken  as  a  perfect  sphere,  and,  for  consistency, 


V  $  =  -g  Ik 


(4) 


where  g  is  taken  as  constant. 
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Only  Eq.  (1)  and  Eq. 
point  of  the  development: 


(2),  in  components,  need  be  expanded  at  thla 


If 


*  A  #*“> +  r  A  w 

+  A  +  2pr 


.  B1 


tan  0 


(5) 


8u  m  u  _ 

8t  “  "  r  cos  0  8X 


8H  .  uv  v  8u  „t  8u 

a7  +  — tan0--g^-w^ 


+  2  Q  v  sin  <f>  - 


FTc 08  0  lx  -  2nwc°8  0-  u7  (6) 


8v 

8t 


§w 

at 


u  8v  v  8v  8v 

~  r  cos  0  0X  ~  r  80  ~  w  dr 

-  2  Ou  sin  0  -  fljl:  -  v  - 
pr  80  r 

u  8wr  y  8w  8w 

"  r  cos  0  ax  "  r  80  ~  w  8r 


+  u2  +  v2  _  1  ?E. 

r  ’  p  Br 


Y  tan  0 


2  Ou  cos  0 


(7) 


(8) 


The  requirement  of  symmetry  about  the  polar  axis  implies  that  all 
parameters  must  be  independent  of  the  longitude,  X  .  We  also  make  the 
usual  approximations  which  are  based  on  the  fact  that  the  atmosphere  is 
a  thin  layer  on  the  earth; this  is, incidentally,  the  basis  for  Eq.  (4).  Th- 
terms  in  w/r  ,  which  arise  due  to  the  expansion  of  a  sphere  with  r  ,  are 
neglected.  Where  r  appears  outside  of  the  differentials  it  is' replaced 
by  a  ,  the  radius  of  the  earth. 

The  Coriolis  term  arising  from  w,  i.e.  2&WCOS0,  is  neglected. 

The  hydrostatic  approximation  is  made  by  replacing  Eq.  (8)  with  the 
hydrostatic  relationship. 

The  coordinate  y  =-’  a0  is  introduced  where  <p  is  expressed  in 
radians.  The  coordinate  r  is  replaced  by  z  (height  above  sea  level' 
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in  the  differentials.  Thus,  our  system,  with  the  inclusion  of  Bq.  (3), 
which  is  so  far  undeveloped,  becomes 


If-  -  ^£<',vco,*>  *  ®<‘>w) 


Sf-  -  jsh  !<“'»•*>  -  wf«  +  20v,to* 


8v  8v 

'  "  w 

-  ju2  +2nusln<>+i^-| 


-  -  &  -  K 

p  8z  8 


It  is  our  preference  to  use  the  hydrostatic  relationship,  Eq.  (12), 


in  integrated  form: 


00 

=  Sl 


The  hydrostatic  pressure-tendency  equation  may  be  obtained  by  differ¬ 
entiating  Eq.  (13)  with  respect  to  time,  and  substituting  for  dp/dt 
according  to  Eq.  (9).  This  yields 


= 

8t 


COS  <f> 


00 

/I* 


v  cos  <t>)  +  gpw 


We  have  yet  to  derive  the  pertinent  Richardson's  equation,  which  is 
now  implicit  in  our  system  because  of  the  stringency  of  the  hydrostatic 
relationship.  In  the  hydrostatic  system,  it  falls  on  w  to  maintain 


4* 


this  relationship.  «•  first  expand  Eq.  (3) : 


£  jl?  *  •  i*"}  •  £  j  I?  *  *ff  *  *&} 

We  then  substitute  for  8p/8t  according  to  Eq.  (14),  and  Bp/Bt 
according  to  Eq.  (9)  : 

?  ^  <p*  ^  + 


p  |  cos  <f>  9y 


1  (pv  cos  <p) 


-  £  (pw)  +  wll  +  V-3? 


8z 


az  ’  ay 


M.  Q 

V  Q  * 


The  latter  equation  is  then  developed  into  the  fora 


dw 

BZ 


1 

P 


K  COS  <t> 


so 

f  ^(pvcos  *)  dz  -  ^(pvcos*) 


+  r  v  f£  +  Pc’  • 


(15) 


The  inclusion  of  the  derived  Eq.  (15)  into  our  system  makes  a  one- 
equation  redundancy.  We  choose  to  dismiss  the  thermodynamic  equation. 
This  choice  is  discussed  at  length  in  Ref.  3.  The  selection  in  this 
case  suggests  later  advantages  in  the  design  of  the  numerical  analogue 
of  the  system.  The  question  of  resolution  is  not  involved  in  this 
choice,  since  resolution  has  already  been  determined  by  the  system  as 
a  whole. 

It  should  be  kept  in  mind  that  with  the  accomplishment  of  a  complete 
system  of  equations,  with  boundary  conditions,  and  based  on  approximated 
laws  and  more  empirical  relationships,  the  mathematical-physical  system 
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xxsil  as  j«(«  • 


has  been  wholly  defined.  Transformations  of  the  equation  system, 
including  changes  of  dependent  and  independent  variables  are  in  them¬ 
selves  completely  without  consequence  and  are  rather  a  means  to  such 
advantages  as: 

(1)  Pons  which  are  more  revealing  for  understanding  of  the 
physical  processes 

(2)  Forms  which  are  more  suitable  for  the  introduction  of 
modifications 

(3)  Forms  which  are  more  revealing  as  to  functional  approach, 
including  indigenous  functions  and  power  series  solutions 

(4)  Forms  which  are  suggestive  in  leading  more  directly  to 
better  numerical  analogues  in  terms  of  ease  of  calculation, 
error  control,  and  the  nature  of  the  analogue. 

At  this  point  we  are  exhibiting  our  system  by  Eqs.  (9),  (10),  OD, 
(13),  and  (15).  One  more  transformation  is  now  carried  out  which  will 
take’ advantage  of  the  special  features  of  the  axially  symmetric  mechanics 

Introduce, 


co 


=  a  + 


u 

a  cos  <p 


(16) 


T)  —  p  COS  <f> 


(17) 


v 


p  CO  cos  <p 


(18) 


where  co  is  the  absolute  angular  speed,  and  2tt  a*)  and  2ir  a  v  are, 
respectively,  the  mass  and  absolute  angular  momentum  of  a  zonal  ring  of 
unit  meridional-plane  area.  With  these  parameters,  Eq.  (9)  transforms 

directly  into 


_a_ 

at 


v 


(19) 
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Equation  (10),  combined  with  Eq.  (16)  leads  to 

•ft  v  "  -  (v*')  -  5F  «0) 

The  parameter  w  replaces  u  in  Eq.  (11) : 

dv  8v  8v  1  8p 

8T  v  8y  ‘  w  8 Z  -  Pfy 

2  2 

-  a  cos  <p  sin  4>  (t*>  -  0  )  .  (21) 

These  equations,  Eqs.  (19),  (20),  and  (21),  together  with  Eqs.  (13)  and 
(15\  are  an  expression  of  the  systea. 

It  should  also  be  stated  that  w  *  0  at  z  =  0  (smooth  earth)  and 
v  =  0  at  the  pole,  where  <P  »  r/2. 

In  these  experiments  we  have  extended  the  model  over  only  one 
hemisphere,  and  imposed  the  artificial  condition  that  v  ■  0  at  <p  =  0  . 
In  effect  we  have  a  wall  at  the  equator.  The  upper  boundary  conditions 
require  that  the  momentum  remain  bounded.  The  effective  boundary  con¬ 
ditions  are  revealed  by  the  discretization,  and  by  the  numerical  analogue. 
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III  THE  NUMERICAL  ANALOGUE 

The  hemisphere,  as  represented  by  distributions  in  a  representative 
meridional  plane,  is  divided  by  surfaces  y  *  ym  ,  and  z  *  *n  , 
which  create  zonal  rings  as  volume  modules.  The  m,nth  volume  module 
consists,  in  the  meridional  plane,  of  one  corner,  two  sides,  and  its 
interior.  All  parameters  referred  to  the  m,n  module  are  subscripted  by 
m,n.  The  principal  evolving  parameters  are  located  at  points  shown  in 
Fig.  1. 


FIG.  1  THE  ARBITRARY  m.n  VOLUME  MODULE 

i 


The  positioning  of  the  parameters  comes  about  in  the  design  of  the 
numerical  analogue.  In  this  way,  the  fields  are  approximated  by  vectors. 
We  shall  circle  the  components  of  the  vectors  in  order  to  differentiate 
them  from  interim  uses  of  the  parameters  in  the  development  of  the 
numerical  analogue. 

The  ym  increase  in  regular  intervals,  6y  =  a7r/98  ,  from 

yQ  =  0  ,  at  the  equator,  to  y^g  =  B.ir/2  ,  at  the  pole.  Levels  are 

selected  in  irregular  intervals,  from  zq  =  0  at  sea  level  to 

z1Q  «  37  km  at  the  base  of  the  uppermost  layer  of  volume  modules. 

18 

The  levels  were  selected  according  to  some  rather  limited  statistics  on 
variability  in  the  atmosphere.  These  statistics  are  not  particularly 
pertinent  to  the  present  experiment. 
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For  the  formulation  of  upper  boundary  conditions,  it  will  be  seen 
that  there  is  little  difference  between  extending  the  upper  layer  of 
modules  from  z  »  37  km  to  z  *  00 ,  or  terminating  the  layer  at,  say, 
z  *  60  km.  The  difference,  which  cones  into  modeling  the  module 
dimensional  constants,  lies  within  the  truncation  error,  and  does  not 
affect  the  degrees  of  freedom.  Ve  have  labeled  the  top  as  z^g  -  60  km  . 

Hie  arbitrary  m,n  volume  module,  as  depicted  in  Fig.  1,  also 
contains  the  lines  y  -  ym_j./2  an<*  z  =  zn*  *  ^ereas  ym_i/2  lies 
midway  between  ym  and  ym_j  >  some  statistics  are  brought  in  in 
selecting  The  statistics  come  in  the  form  of  the  N.A.C.A. 

Standard  Atmosphere.  The  Standard  Atmosphere  values  of  the  state 
parameters  are  denoted  by  the  subscript  s  .  The  level,  zq4,  ,  is 
determined  by  the  relationship 


i 

J 


n 


psdz  *  psn*  (zn  -  zn-l> 


(22) 


n-1 


The  level  at  which  Pg  *  Psn*  is  the  level  zn*  ,  where  the  standard 
density  is  a  mean  for  the  layer. 

With  the  discretization  of  the  field  established,  we  proceed  with 
the  design  of  the  numerical  analogue  beginning  with  the  simplest  of  our 
equations,  the  hydrostatic  relationship,  Eq.  (13).  The  pressure  i.3 
simply  the  weight  of  the  column  of  air  above.  This  is  developed  as 
follows: 


19 


®m,n  ’  *  /  p  dz 


n’ 


=  g 


I 


n  19 

p  dz  + 


19  r 

l  I 


p  dz 


n1 


r=n+l  z 


r-1 
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Here  r  is  used  as  a  dummy  variable  for  n  .  Thus,  we  aake  the 
approximation 


‘  \  2  ®m.n  +  I  ™r®m,r 

[  r^i+1 

in  which  our  firat  module  dimensional  constant  is 


z 

n- 


1 


(23) 


(24) 


In  order  to  conserve  symbols,  all  module  constants  will  be  denoted  by 
boxed  numbers,  subscripted  according  to  their  dependency.  Further 
statistical  refinement  of  the  module  constants  is  not  warranted  at  this 
time,  but  should  be  considered  in  later  developments  and  applications 
of  the  model. 

Equations  (19)  and  (20)  are  both  of  the  form 

-L  *  =  V  .  <  *  V  )  (25) 


in  a  cartesian  y,z  frame,  where  *  =  rj  and  v  ,  respectively. 
Analogues  for  these  are  ideally  constructed  by  first  integrating  over 
the  y,z  area  of  the  volume  module,  and  applying  Gauss's  theorem: 

-|r/*da  =  -/*V'df  '  (26) 

m  A  L 

That  is,  the  total  change  rate  in  the  module  is  equal  to  the  net  flux 
rate  into  the  module.  In  designing  the  analogues,  one-to-one  correspond¬ 
ences  are  imposed  between  the  mass  of  the  module  and  the  density  at  an 
interior  point,  and  between  the  absolute  angular  momentum  of  the  module 
and  pu>  at  an  interior  point. 


The  aas a  of  a  voIum  module,  omitting  the  factor  2*a  ,  ia 
approxiMtad  by 

,?n  ym 

J  j  pcoa^dyto  *@mnmnQIm  <27) 

*n-l  ^m-1 

and  tha  absolute  angular  aomentua  of  a  volume  module,  omitting  the 
factor  2«a2  ,  is  approxiMtad  by 

poi  cos3  0  dydz  *  (pw)^  (28) 

Vl  ym-l 


where 


J  cos  0  dy 
ym-l 


m 


QJm  =  f  cos  0  dy 
m 

ym-l 


(29) 


(30) 


For  the  fluxes  through  module  interfaces,  the  speeds  are  suitably 
located.  However,  p  and  pu>  must  be  interpolated.  Linear  inter¬ 
polations  are  used  except  for  p  in  the  vertical,  for  which  the 
Standard  Atmosphere  is  again  invoked. 

For  p  ,  at  the  level  ,  is  used 
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where 


-  i. 


*n+l,*  "  *n*  ps.n+l,* 


C5LS 


n  z 


zn+l  *  ~  gn  pan 
n+1,*  -  *n*  psn* 


The  interpolation  ia  such  that  it  would  be  exact  if  the  diatrlbution 
were  Standard  Ataoaphere. 

The  aass  flux  ratea  are  repreaented  by  n  »nd  ®m  n  ,  *nd 

the  abaolute-angular-aoaentua  flux  ratea  by  (c)  and  (5}  ,  aa 

m,n  v-^m,n 

ahown  in  Fig.  2. 


y  z 

Jm  n 


®  .  ® 


FIG.  2  THE  FLUXES 


In  this  manner,  that  is  by  first  integrating  over  the  volume  module, 
Eqs.  (19)  and  (20)  are  approximated  by 


£®mn  =  <»«> 

i  ©nn  *  ffln®m{©m-l,»-©™,+®m,»-l-®J  <35> 


WK-VM...-. 


where  QOn  □□  m»nd  CO  m  are  the  inverses  of  CDnCZ!)m  and  C0m  re¬ 
spectively.  The  fluxes  are  calculated  according  to 

®mn  *  ©ran  j®m+l.n  *  ®mnj  * 

©ran  S  ©ran  {  ®m.n*l  *  ™m  <”> 

©mn  5©mn  |©m+l,n  +  ©  mnl  “»  00 m  ' 


where 


\  m  rm 

mn/  n  m 


®mn  ^®n  ®m(lH-l  +  ®mn| 

{™„  ©m,»l  *  ™n  ©mn 


rm  m  (39) 


2  cos 


™  m  *  co"  *m-l/2  6y 


rm  =  ^  cos3  <pm 
u"-*  m  2  m 


n  nT 
zn+l,*  "  zn* 


zn+l,*  ~  Zn 
zn+l,*  zn* 


™m  s  cos  *m-l/2  6y 


wswsbwbmbss* 


All  fluxes  mre  through  interior  Interfaces.  They  transport  froa 
one  aodule  into  the  other.  It  is  apparent  that,  because  the  healsphere 
is  contained,  the  numerical  analogue  strictly  conserves  (except  for 
random  round-off)  total  mass  and  total  absolute  angular  momentum: 

£  CD  (p)  *  constant  (46) 

‘"n  m  vC/mn 

1  m  m  (pH?) «  constant.  (47) 

n  m  v /mn 

The  suas  are  taken  over  all  modules.  These  sums  are  carried  out  for  each 
tlae  step  during  the  Integration  as  computational  checks.  It  is  also 
noteworthy  that  vertical  mass  fluxes,  alone,  do  not  change  sea-level 
pressure. 

The  analogue  for  the  Richardson's  equationlEq.  (15),  is  also 
constructed  by  integration  over  the  elemental  m,n  volume  nodule.  How¬ 
ever,  before  Integrating,  we  wish  to  allow  for,  and  make  explicit,  some 
dependency  of  p  and  p  on  z  ,  in  the  module.  We  approximate  the 
statistics  of  this  dependency  by  isothermal  conditions  at  temperature 


P  =  P„* 


(48) 


P 


(49) 


where 


E 


RT 


<z  •  v> 


sn’ 


(50) 
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+  “n  Um  (®mn  - 
+  Ll2J  |  0m-i,n  +  ©mnj  |(£)m+l,n  ”  ®m-l,n 


+  6y  Q, 


mn 


(54) 


This  defines  some  additional  parameters: 


as: 


g 


m  r  cos  0m.l/2 


(55) 


rnn  =  -  T  „ 

g  sn" 


an  =  -  t  „ 

g  sn' 


'  RT 


—  T  *  \{ 
g  sn*  I  ' 


1  -  e 


fc  |vVl|l 


z  -  z  , 
n  n-1 


z  -z 


1  -  e 


RTsn*  1  n  "-1 


}  (56) 


(57) 
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r 


ran 


m  cos 


m-1/2 


(58) 


V' 


l 

l 

i 


£  I 


nffl  ■  R/4C 


mn 


R 

C  ^sn*  ^nin 
P 


(59) 


(60) 


In  the  forcing  torn,  Bq.  (60),  we  have  wade  the  approximation  of 
replacing  ®m>n  by 

psa*  ’ 

The  remining  equation  for  which  an  analogue  ia  desired  ia 
Bq.  (21) .  The  analogue  ia  constructed  by 

df  ©inn  ”  "  ©mn|©n+l,n  ~  ©m-l,n  J  2  6y 

-  a^SU4®mtl,/  ®m-l,n-l  *  ®m,n-l 


[/ — \  /~\  \  ®m+l,n  ~ 

(©m,n*l  - 


©m+l,n  '®mn  2_ 

6y 

m+l,n 


-  eh 


m 


(®m*l,n  *  © 


mn 


2-4Q2 


(61) 


The  additional  parameters  are 


EEL  = 


n  4<Vl  *  -  Vl  y 


EHm  = 

m 


a  cos  sin 


(62) 


(63) 
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The  space  dependency  of  our  mathematical-physical  system  is  thus 
numerically  modeled  by  Eqs.  (23),  (34),  (35),  (54),  and  (61).  Also 
implicit  in  our  system  is  the  perfect  gas  law,  which  yields 

©nui  ‘  ®m«  /  R®mn  0 

Some  further  discussion  of  the  boundary  conditions  is  necessary. 

The  lower  boundary  condition  is 


(65) 


The  vertical  motion  is  thus  diagnosed  by  progressing  upward  with  the 
analogue  of  Richardson's  equation,.  Eq.  (54). 

The  upper  boundary  conditions  essentially  take  the  form 


®m,19*  ®  ■».!!.*  0 


(66) 


This  implies  that  we  consider  z to  be  effectively  at  infinite 

height.  Examination  reveals  that  the  heating  in  the  upper  layer  of 

modules  has  little  effect  on  the  evolution  via  the  analogue  of 

Richardson's  equation,  Eq.  (54).  Because  of  Conditions  (66),  the 

calculated  (w)  are  used  only  in  Eq.  (61)  ,  in  determining  a 

— '  m ,  iy 

mean  w  for  the  vertical  advection.  The  influence  there  is  rather 
minor . 

More  realistic  influence  of  the  upper  layer  heating  may  be  built 
into  the  model  as  follows.  Equation  (23)  yields 

®m,  19  =  g  ^  19(2  ©m,!^  ( 
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This,  in  effect,  specifies  the  temperature  there.  According  to 
Eq.  (65),  then. 


®m,19  “  flR  ®  19  ’ 


(68) 


By  preacribing  03  19  as  a  function  of  m  and  time,  the  heating  of 
the  upper  layer  nay  be  modeled  into  the  analogue.  This  requires 
suitable  modeling  considerations.  In  our  experiments,  CH 19  »• 

left  as  a  constant. 

The  boundary  conditions  at  the  pole  and  at  the  equator  are, 
respectively, 


(69) 


(70) 


Additional  boundary  specifications  are  required  because  of  the  use  of 
some  double  differences  in  Eqs.  (54)  and  (61).  In  Eq.  (54),  in  the 
term  with  the  coefficient  CEI  ,  appears  the  centered  double-interval 

difference: 


®m+l,n  ”  ®m-l,n 


(71) 


For  m=s49  this  is  replaced  by  the  one-sided  difference 


(72) 


and  for  m 


1  It  is  replaced  by  the  one-aided  difference 


*(©*.-©...)  •  <ra> 
Xn  Eq.  (61),  the  centered  double-interval  difference 

^Dm,n+1  ©m,n-l 

is  similarly  modif ied  for  n  ■  1  and  n  ■  19  .  This  makes  the  numerical 
analogue  system  a  complete  system. 

We  have  followed  the  procedure  of  first  numerically  modeling  the 
space  dependency  of  our  system.  It  remains  to  discretize  the  time 
dependency  and  to  indicate  the  integrating  procedure.  The  time  con¬ 
tinuum,  t  ,  is  discretized  by  time-increment  fit  ,  such  that 

t  *  t  fit  ,  T  »  0, 1,2,3  . ..  <74> 


where  T  is  the  number  of  the  time-step. 

Initially,  and  at  any  subsequent  time-step,  the  complete  specifi- 


cation  shall  consist 

of 

the  principal 

dependent  parameters 

\ 

0mn 

for 

m 

*  1,2.  ..49  ; 

n  =  1, 2. . .  19 

®mn 

for 

m 

*  1,2.  ..49  ; 

n  *  1,2...  19  ) 

©  mn 

for 

m 

=  1,2.  ..48  ; 

n  *  1, 2. . .  19 

From  these  we 

may  derive 

(76) 
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and  diagnose^m  q  according  to  Eq.  (23)  and(w)m  Q  according  to  Eq.  (54). 

The  principal  dependent  parameters,  (75),  at  t  *  r  fit  ,  may  be 
represented  by  the  vector  V  .  Equations  (34),  (35),  and  (61)  give  us 


i.  v 

dt  ~ 


as  non-linear  functions  of  the  parameters.  This  derivative  is  approxi¬ 


mated  by  the  double-interval  difference 


|Vl  -  Xt-i  )  /2  «* 


which  gives  the  conventional  "leap-frog"  integration  procedure. 


V  ..  =  V„  ,  +  2  fit  V 

~T+1  ~T-1  dt  •'* 


This  completes  the  numerical  analogue.  To  begin  the  time  integration, 
using  Eq.  (79),  requires  that  we  know  for  two  consecutive  time- 

steps. 

Stationary  solutions  of  the  numerical  analogue,  i.e.,-^  V  =  0  , 
include  the  class  of  zonal  flows 


with  ®m,n  and  Gln.n  satisfyinS 

©m+l,n  ~®mn  2 
®mn  +  ®m+l,n  ^ 


+  ‘21)m  (®m+l,»  +©m„)2-  41,21  " 


=  0  . 


Equation  (81)  derives  f roe  Bq.  (61) .  Equations  (80)  and  (81)  also 
include  the  degenerate  case  of  relative  rest: 

©mn  '  0  •  ©mn  E  0  •  ®mn  “  “n  ' 

In  these  experiments,  the  integration  is  begun  by  specifying  a 
stationary  solution,  representing  conditions  both  at  T  0  and 
T  «  1  .  The  equilibrium  is  upset  by  the  commencement  of  heating  at 
T  «  0.  The  heating  is  introduced  in  steps,  as  discussed  in  Sec.  I. 

The  experiments  discussed  herein  are  limited  to  integrations 
begun  with  a  Standard  Atmosphere  at  relative  rest.  Heating-rate 
distributions,  Eq.  (60), 

q  for  m  *  1  2.  ..49  »  n  *  1,2.  ..19 
Mmn 

are  specified  constant  in  time,  except  that  they  are  gradually  intro 
duced  as  follows: 


Time 

Interval 

Proportional  Rate 

0 

to 

2  at 

1/32 

6t 

to 

3  at 

1/16 

2  at 

to 

4  at 

1/8 

3  at 

to 

5  at 

1/4 

4  at 

to 

6  at 

1/2 

5  at 

to 

7  at 

1 

etc. 

1 
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The  levels  selected  and  the  Standard  Ataosphere  values  for  calcu¬ 
lation  of  the  constants  are  shown  in  Fig.  3.  The  values  and  (t)d 

correspond  to  the  initial  specification  ,  and  are  determined 

by  Eqs.  (23)  and  (65)  respectively.  Additional  values  required  are: 

60  s  0.03205706  radians 
a  =  6,371  k* 

-2 

g  -  9.80665  neters  secs 
R  =  287.04  kj  ton”1  deg  1 
Cp  =  1004  kj  ton  1  deg  1 

The  integrations  are  based  on  the  meters-tons-seconds  systea  of  units. 

The  computational  stability  of  the  analogue  and  the  magnitude  of 
the  time  increment,  fit  f  are  discussed  in  Sec.  IV. 
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XV  TEST  INTEGRATIONS,  COMPUTATIONAL  STABILITY  AND  GRAVITY  WAVES 


The  "leap-frog"  method  of  time  integration,  Eq.  (79),  has  several 
recognized  undesirable  properties  which  limit  its  usefulness.  It  is 
used  here  in  the  interests  of  simplicity  and  economy. 

The  centered  single-difference  analogue, 


d_  1 
dt  2 


(83) 


is  more  accurate  and  apparently  more  natural.  It  is  implicit  and  may 
be  used  by  reiterating  each  time  step.  The  convergence  condition  imposes 
an  upper  limit  on  fit. .  The  programming  is  more  complex  and  the  compu¬ 
tations  are  lengthier  than  for  the  "leap-frog"  method.  For  these  reasons 
it  was  not  used  in  our  limited  experiments. 

These  analogues,  Eqs.  (79)  and  (83),  are  readily  analyzed  for 
(7) 

linear  systems'", 


_d 

dt 


Y  -  *Y 


(84) 


where  the  matrix  M  is  an  analogue  of  the  linear  operations.  Such 
analysis  is  quite  pertinent  to  non-linear  systems  because  of  the  high 
degree  of  inherent  linearity  exhibited  by  evolutions  and  phenomena 
which  are  governed  by  such  systems.  This  is  especially  so  for  gravity 
waves  in  the  atmosphere,  as  will  be  demonstrated  later. 

Our  numerical  analogue  was  programmed  for  a  Burroughs  220  computer 
having  a  5000-word  core  memory  and  adequate  magnetic  tapes.  The  magni¬ 
tude  of  our  problem  required  a  considerable  amount  of  data  transfer 
between  tapes  and  cores  for  each  time-step.  As  programmed,  the  computer 
running  time  per  analogue  time-step  required  approximately  1.25  minutes 
of  calculations  and  1.15  minutes  for  tape  movement.  In  addition  the 
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full  output  of  the  field*  ®  ,©  »  ©  »  ©  »  fe  »  and  ®  for  a  8poci" 
fled  time-step  required  approximately  2.75  minutes  of  on-line  tabulation. 
With  output,  thi*  total*  to  an  approximation  of  5.15  minutes  per  tiae- 

•tep. 

The  *ua  check,  for  total  mass,  Eq.  (46),  and  total  absolute  angular 
aoaentua,  Bq.  (47),  are  calculated  for  each  tiae  atep. 

The  heating  rate,  Q  ,  ha*  the  diaen.ion*  lenergy  aaa*  1  tiae  1 1  , 
whether  considered  as  parcel  or  local  input.  In  our  *y*tea  the  Q  field 
is  locally  specified. 

A  frequently  encountered  convention  is  to  express  heating  in  terms 
of  the  equivalent  teuperature  change  of  a  parcel  at  constant  pressure: 


b  - 


cpte 


(85) 


The  interpretation  of  this  specification  can  be  misleading  if  it  suggests 
its  tiae  integration, 


6Q  *  C  6T 

F 


(86) 


which  is  not  generally  meaningful  in  the  circulating  atmosphere.  The 
full  relationship  for  local  temperature-change  rate,  implicit  in  our 

system,  is 


Q 


v-^T-^ar-p  v  ^ 


(87) 


In  the  case  of  a  column  whose  motion  is  restricted  to  thermal 
expansion  in  the  vertical,  each  fluid  layer  remains  at  constant  pressure 
hydrostatically.  For  this  evolution,  Eq.  (87)  reduces  to 


Q 


=  cP!§r 


9T 
+  w  Wz. 


(88) 


32 


4. 


and  Bq.  (15)  reduces  to 


This  experiment  was  designed  in  order  to  excite  the  •■elleet-ecale  modes, 
with  the  prospect  of  physically  interpreting  the  evolution. 

In  order  to  anticipate  the  evolution,  we  ahould  expect  the  following 
sequence  of  events:  The  heated  column  expands  upwards.  This  causes  an 
excess  of  pressure,  over  adjacent  columns,  at  all  levels  above  the  sur¬ 
face.  The  horizontal  pressure  gradients  so  established  accelerate  the 
air  out  of  the  heated  column  both  northward  and  southward.  Because  of 
this  upper  mass  divergence,  a  pressure  deficit  first  appears  at  the 
bottom  of  the  heated  column  and  is  progressively  shown  by  a  deeper 
surface-based  layer.  This  causes  a  lower  mass  convergence  and  the  cir¬ 
culation  pattern  becomes  established. 

For  trial  run  number  one,  fit  =«  1  hour  was  arbitrarily  selected. 
Computational  instability  was  expected  for  so  large  a  time  increment, 
and  the  integration  was  carried  out  only  to  time  step  r  =  6  .  The 
computational  instability  was  violent,  because  the  tendencies  were 
allowed  to  act  for  2fit  *  2  hours  without  being  checked  by  their  conse¬ 
quences.  The  horizontal  mass  divergence  appeared  at  upper  levels  in  the 
heated  column,  at  T  *  3  .  By  allowing  it  to  act  from  T  -  2  to 
T  =4  ,  it  more  than  evacuated  the  upper  three  modules,  producing 
negative  density. 

For  trial  run  number  two,  the  time  increment  was  reduced  to 
fit  =  6  minutes,  and  the  integration  was  taken  out  to  x  =  16  •  Again 
we  had  computational  instability,  but  the  growth  rate  was  slow  enough 
to  permit  a  careful  analysis. 

The  computational  instability  was  manifested  as  an  oscillation  of 
mass  between  even-  and  odd-numbered  whole  columns.  The  period  of  this 
oscillation  was  4  time  steps.  The  disturbance  propagated  northward 
and  southward  from  the  heated  column  at  the  rate  of  one  column  every 
two  time  steps. 


For  the  disturbance  growth-and-propagation  mechanism,  only  a  few 
terms  acted  in  the  numerical  analogue.  In  these,  where (p)m+i>n  + 

® m  n  occurred>  it:  could  be  approximated  by  2  Pgn*  .  In  effect  the 
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analogue  reduced  to 


S’  ®mn  " 

(94) 

®mn  * 

©ron  2 p»n* 

(95) 

!•©««,  * 

-1  .  I/C\  \ 

(96) 

Pan+  6y  |v£/m+l,n  viyuml 

and  the  hydrostatic  pressure,  Eq.  (23) .  These  four  equations 

form  a 

complete  linear  aystea. 


Apparently,  neutral  compressioa-gravity-modc  oscillations  between 
adjacent  columns  had  been  transformed  by  the  numerical  analogue  into 
growing  modes.  Since  it  appeared  that  the  disturbance  mechanism  was 
essentially  linear,  the  effect  of  the  time-integration  procedure  was 
linearly  analyzed. 

We  analyze  the  effect  of  approximating 

V  -  M  V  (97) 

Cat  »v  ®  ^ 


by 

V  ,  =  V  ,  +  2  fit  M  V  .  (98) 

~T+1  ~T-1  «  ~T 

The  system  of  Eq.  (97)  has  characteristic  solutions  (modes) 

V(t)  =  Re:  V  e1Ct  =  Re:  V  (e1C6t')T  •  (99) 

Substitution  in  Eq.  (97)  yields  the  system  of  equations 

[  M  -  ic  I  ]  V  =  0  (100) 
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where  |  ia  the  identity  Matrix.  This  imposes 


is  -  *»ii 


(101) 


This  characteristic  polynomial  has  in  general  as  many  roots,  C  ,  as 
y  has  components.  Corresponding  to  each  characteristic  value,  Cj  ,  the 
system  of  equations  in  Eq.  (100)  yields  a  characteristic  vector,  V.  • 

The  system  of  linear  equations  in  Eq.  (98)  has  characteristic 
solutions 


V,  -  Re:  V/ 

iw  9 


(102) 


Substitution  in  Eq.  (98)  yields 


M  -  y2~  1 

«  2  fity 


(103) 


In  comparing  Eq.  (103)  with  Eq.  (100),  we  see  that  they  have  the  same 
characteristic  vectors,  Vj  ,  and  that  corresponding  to  the  characteristic 
value,  C.  f  we  have 


2  fit  y. 
J 


‘  cs 


(104) 


This  quadratic  has  the  two  roots 

y.  =  i  fit  C .  ±  ~\Jl  -  (fit  C.)2  .  (105) 

The  positive  root  approximates  the  proper  modes,  while  the  negative  root 
yields  extraneous  modes  admitted  by  the  "leap-frog"  method. 
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Consider  a  neutral  node  with  characteristic  value  Cj  (!•*•>  Cj 
real)  and  a  tiae  increment  such  that 

6t  Cj  >  1  .  (106) 


The  proper  node  of  the  "leap-frog"  method  yields  the  characteristic 
value 

Vj  -  i  jat  ci  +  s  <107> 

where  p  is  real  and  greater  than  one.  Thus,  the  neutral  node  is  trans¬ 
formed  into  an  amplifying  mode  which  satisfies 


v 

~j.T+2 


This  explains  the  evidenced  growth  cycle  of  four  time  steps. 
All  neutral  modes  for  which 


(108) 


>  l/6t  (109) 

are  transformed  into  amplifying  modes,  and  the  greater  Cj  ,  the  greater 
is  the  resulting  growth  rate.  This  is  what  we  assumed  has  happened  in 
our  integration. 

Equation  (108)  provides  us  with  a  test  which  we  can  apply  to  trial 
integration  number  two.  For  the  component  V  =(^)25,19  ^  turns  out 

that 


V  /V 
v9/  v7 

=  -3. 328 

Vl</V8 

=  -3.329 

VU/V9 

=  -3.493 
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VV10  -  -J  «8 

V13/Vll  '  *S-512 
V14/V12  *  -3'504 

VV13  '  -»-482 

V16/V14  ‘  -3'228 

The  component  V  “©25, 15  ***  al*°  trl,<,: 

V15/V13  -  -3'626 
V16/VM  -  -3-616 

These  ratios  are  reaarkably  steady. 

Since  the  oscillations  occur  between  adjacent  whole  columns,  and 
we  have  49  columns,  there  must  be  a  whole  family  of  modes  with  oscil¬ 
lation  periods  spread  over  a  narrow  range.  This  spread  should  arise 
because  of  the  m-dependency  of  the  module  constants.  The  ratios  above 
may  reflect  mode  predominance  because  of  the  manner  in  which  the  motion 
was  initially  excited. 

As  a  reasonable  upper  limit  on  the  ratios  we  select 


2 


*  4 


This  yields 


fit  C. 


1.25 


and 


C.  «  0.21 
J 


radian  per  minute  . 


(110) 
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This  corresponds  to  the  period  of  oscillation 


7 r-  »  1810  seconds  .  (HD 

Cj 

Another  Method  of  approxiMately  determining  this  characteristic 
period  is  by  the  following  somewhat  crude  analysis  of  a  simplified  model. 
Consider  the  compression-gravity  oscillation  between  two  adjacent  volume 
modules.  Equations  (94),  (95),  (96),  and  (23)  yield: 

i®2  *-&©!  -  in.ni,® 

®  =  mnE)m  2f>sn*© 

S’  ©  “  '(®2  ‘  ®l) /psn*  ^ 

®1  =  *  2  ®1 

@2  g  ^  n  2 

From  these  equations  we  derive 

i  ©  ■  -k2  ® 

where 

k2  =  If  mf  ej„  unm 

k  «  \rg6z  /  5y  .  (112) 
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Thii  further  shews  thst  the  spread  in  periods  which  this  family  of  modes 
has  in  our  model  due  to  a-dependency  is  probably  rather  slight.  As  an 
upper-bound  estimate  we  put  6z  *  25  km  .  This  yields 

k  *  0. 15  radian  per  minute 

which  is  somewhat  lower  than  our  over-all  estimated  upper  bound, 

Bq.  (110). 

The  preceding  analysis  implies  that  this  form  of  computational 
instability  should  be  eliminated  by  choosing  a  time  increment 

fit  <  l/C 

max. 


We  have  estimated  ^max  at  0.21  radian  per  minute.  Thus,  the  upper 
limit  on  at  is 


»  4.8  minutes 

In  all  subsequent  integrations  we  use 

6t  *  3  minutes 

This  means  that,  as  programmed  for  the  Burroughs  220,  the  integration 
progresses  at  about  real  time. 

At  this  point  in  the  experiments  the  computer  program  was  speci¬ 
fied  to  tabulate  a  set  of  outputs  at  intervals  of  20  time  steps.  Such 
a  set  of  outputs  consists  of  complete  tabulations  of  each  of  the  six 
fields,  P,T,w,p,o;,  and  v  .  The  purpose  in  having  two 
successive  values  in  time  is  to  inspect  for  stability  and  smoothness 
of  the  evolution,  in  particular,  to  determine  the  growth  of  extraneous 
behavior  admitted  by  the  "leap-frog"  method.  The  availability  of  the 
fields  in  high-speed  memory  favored  the  output: 
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Fields  P,  T,w,  for  r  «  1,2;  21,  22;  41,  42;  etc. 
Fields  p,i i),v,  for  r  *  2,3;  22,  23;  42,  43;  etc. 


(113) 


Thus,  we  have  all  six  fields  for  t  *  2,22,42,62,  etc.  In  addition,  the 
sua  checks  for  mass  and  absolute  angular  aoaentua  were  tabulated  for 
every  tiae  step. 

It  should  also  be  mentioned  that  the  coaputations  are  carried  out 
in  "floating  point"  arithmetic  with  approxiaately  eight  decimal  digits 
carried  throughout. 


V  RESULTS,  ONE  COLUMN  OF  HEATING 

A  Standard  Atmosphere  at  relative  rest  is  disturbed  by  tbe  steady 
heating  field  prescribed  by  Eqs.  (93).  This  heating  is  Halted  to  the 
region  between  two  latitudes  (Column  25  in  the  representative  meriodlonal 
plane)  at  mid-latitude.  The  equivalent  temperature  rate,  Eq.  (85),  of 
the  heating  is  uniform  (10°C/24  hours)  at  all  levels  in  the  region. 

The  heating  is  commenced  in  steps,  as  discussed  in  earlier  sections. 

following  the  trial  runs  that  were  analyzed  in  Sec.  IV,  the  experi¬ 
ment  was  repeated  using  the  time  increment  of  3  minutes.  The  integration 
was  carried  out  to  t  =  23  (approximately  two  periods  of  the  shortest 
mode  in  our  model).  No  computational  instability  was  discernible. 

Although  this  experiment  was  designed  primarily  as  a  test  of  the 
model,  the  results  appear  to  be  sufficiently  interesting  for  inclusion 
in  this  report.  The  circulation  established  at  T  =22  is  shown  in 
Figs.  4  to  10.  These  figures  show  that  section  of  the  representative 
meridional  plane  over  which  the  disturbance  has  spread  appreciably.  The 
vertical  height  scale  is  linear  in  terms  of  Standard  Pressure,  except 
for  the  upper  two  module  layers  where  this  scale  is  magnified  by  a 
factor  of  ten  (see  Fig.  3).  Because  of  this  change  in  scale,  the  isoline 
analyses  are  not  carried  to  these  upper  layers.  The  equator  is  at  yQ 
and  the  pole  at  y49  ' 

Attention  is  drawn  to  the  fact  that  Figs.  5,  6,  8,  and  9  show 
anomaly,  or  change  from  the  initial  values.  This  is  indicated  by  the 
subscript  "a."  In  the  case  of  Fig.  8, 


U)  =  U) 

a 


-  u/a  cos  <t> 


(114) 


These  figures  speak  well  enough  for  themselves.  Of  particular 
interest  is  the  somewhat  surprising  phenomenon  of  local  cooling  in  the 
one  column  which  is  receiving  thermal  energy.  This  is  due  to  the 
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FIG.  8  ONE-COLUMN  HEATING,  u,a  of  r  *  22 
Units:  10"’  0  rodians  sec"' 
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FIG.  9  ONE-COLUMN  HEATING,  T0  ot  r  =  22 
Units:  10'2  °C 
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established  circulation  which  brings  aass  into  the  coluan  below  and  out 
of  the  coluan  above.  Hie  associated  ascent  and  adiabatic  cooling 
apparently  is  sufficient  to  overcosw  the  heating,  in  local  influence. 

Figure  10  is  a  streamline  schematic.  The  vectors  show  the 
streamline  slope,  with  an  exaggeration  of  100  in  the  w  to  v  ratio. 
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VI  RESULTS,  HEMISPHERE  WINTER  HEATING 

For  our  ujor  experiment,  a  heating  rate,  designed  to  approximate 
mean  winter  conditions,  was  used  as  the  forcing  field. 

For  the  net  radiative  contribution  of  the  heating  we  used  the 
distribution  of  Coulson^1^.  To  this  was  added  an  estimate  of  the  contri¬ 
bution  of  latent  and  sensible  heat  carried  into  the  atmosphere  from  the 
surface  layer,  Fig.  11.  For  the  latter,  an  equivalent  temperature  rate 
was  calculated  on  the  basis  that  this  energy  goes  into  uniform  heating 
of  Coulson 's(l)  model  of  the  troposphere  at  each  latitude.  The  total 
heating  field  in  terms  of  equivalent  temperature  rate,  Tg  ,  is  given 
in  Fig.  12.  This  distribution  was  used  as  the  forcing  field. 


FIG.  11  ESTIMATE  OF  LATENT  AND  SENSIBLE  ENERGY  RELEASE 
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FIG.  12  APPROXIMATION  OF  MEAN  WINTER  HEATING 
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This  approximation  of  asan  winter  hasting  was  introducad  in  stops 
to  disturb  a  Standard  Atnosphora  at  ralativa  rest.  A  tins  incraaont  of 
3  ainutas  was  usad.  Tha  intagration  was  carriad  to  T  ■  183  Tha 
output  was  as  spaciflsd  by  Bq.  (113) . 

Tha  resulting  evolution  was  smooth,  as  shown  by  selected  v  coa- 
pononts  given  in  Table  I  and  graphed  in  Fig.  13.  Careful  analysis  of 
successive  values  in  Table  1  shows  that  very  little  extraneous  behavior 
has  been  admitted  by  the  "leap-frog"  aethod. 

Figures  14  to  19  show  the  distributions  of  w  ,  p  ft  ,  pft  ,  V  , 

u  and  T  respectively,  at  time  step  r  =  22  (approxiaately 
A  A 

Hour  1)  .  The  subscript  "a"  indicates  anomaly,  or  change  from  the  initial 
values. 

The  reaaining  figures,  Figs.  20  to  27,  refer  to  r  =  182  (approxi¬ 
mately  Hour  9) .  Figure  25  is  obtained  from  Fig.  24  according  to 
Eq.  (114) 


u  *  a  cos  0  h>a 


(115) 


Figure  27  is  a  streamline  schematic.  The  vectors  show  the  streaaline 
slope,  with  an  exaggeration  of  100  in  the  w  to  v  ratio. 

The  evolution  is  quite  amenable  to  physical  reasoning.  The 
atmosphere  shrinks  by  cooling,  which  is  stronger  northward.  This 
shrinking  results  in  a  northward  pressure  gradient  and  northward 
motion  at  higher  levels.  As  this  flow  progresses,  carrying  mass 
northward,  a  southward  pressure  gradient  and  southward  motion  develop 
at  lower  levels.  A  major  single  cell  meridional  motion  develops  with 
minor  superimposed  cells.  By  Hour  9,  a  core  of  westerly  current  has 
developed  at  about  50°  North  and  a  height  of  about  10  km. 

The  sum  checks  for  mass  and  absolute  angular  momentum  remain 
steady  except  for  minor  round-off  which  appears  to  be  quite  random. 

As  shown  by  Fig.  13,  the  evolution  is  in  a  stage  of  continuing 
development  at  Hour  9.  The  main  cell  seems  to  be  shifting  and 
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FIG.  13  WINTER  HEATING,  THE  EVOLUTION 
OF  SELECTED  v  COMPONENTS 
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FIG.  14  WINTER  HEATING,  w  ot  r  •  22  H  hour) 
Units:  10"4  cm  sec"' 
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FIG.  16  WINTER  HEATING,  p0  ot  r  -  22  M  hour) 
Units;  10'3  mb 
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FIG.  19  WINTER  HEATING,  T0  of  r  -  22  H  hour) 
Units:  KT3  °C 


FIG.  20  WINTER  HEATING,  w  at  r  =  182  (-9  hours) 
Units:  10'3  cm  sec'1 
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FIG.  24  WINTER  HEATING,  oi  r  .  182  (-9  hours) 
Units:  10'1 0  rodions  sec"' 
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FIG.  26  WINTER  HEATING,  T0  ot  r  ■  182  (-9  hours) 
Units:  10'J  ”C 
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intensifying.  Tbs  predominant  northward  motion  nw  to  reach  *  peak 
at  About  lour  6;  thereafter  tho  southward  flat  begins  to  dominate.  This 
observation  stages t a  that  s  hemisphere  compression-gravity  oecilUtion 
with  period  of  about  14  hours  ia  aapsrlapoaad  on  tbs  ns  rid  ions  1  osllular 
olroulatlon. 

This  major  (rarity  oscillation  oeeurs  bscauso  of  tbs  rslativs 
abruptnssa  of  tbs  winter  beating  applied  to  a  Standard  Atmosphere  at 
relative  rest.  Its  reversals  are  dictated  by  aero  flow  at  the  pole  and 
serosa  the  eguator.  To  avoid  this  oscillation,  the  beating  would  have 
to  be  Introduced  gradually  over  a  period  of  a  day  or  ao. 


VII  RECOMMENDATIONS 
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On  appraising  the  model  and  the  experimental  findings,  a  number  of 
modifications  and  improvements  may  be  put  into  focus. 

For  the  investigation  of  planetary  circulations  driven  by  seasonal 
heating,  the  most  obvious  and  major  improvement  is  to  extend  the  model 
from  pole  to  pole.  This  not  only  removes  the  artificial  boundary 
influence  on  the  hemisphere,  but  also  allows  us  to  investigate  the 
important  seasonal  exchanges  of  mass,  momentum,  and  energy  between  the 
northern  and  southern  hemispheres. 

Another  major  consideration  is  that  the  model  be  reprogrammed  for 
a  faster  computer  with  greater  storage.  This  would  make  experiments 
more  economical  on  the  one  hand,  and  would  allow  the  evolutions  to  be 
carried  to  great  lengths  on  the  other  hand. 

In  planning  to  so  extend  the  evolutions,  additional  physical 
processes  should  be  included.  Those  to  be  considered  are: 

(1)  The  surface  layer  should  be  frictionally  coupled  to  the 
earth  with  a  flow  of  angular  momentum  across  the  lower 
boundary . 

(2)  There  should  be  mixing  in  the  vertical,  attributed  to  con¬ 
vection.  This  may  be  handled  by  exchanges  of  mass,  momentum, 
and  energy  between  modules  along  the  vertical.  This  can  be 
done  periodically  and  should  be  dependent  on  vertical  structure 
and  stability  considerations. 

(3)  A  horizontally  acting  friction  may  also  be  necessary.  A 
friction  of  the  type 

k  (PVH) - +  *  m2vhvh-  (p  V H>]  016) 
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may  be  considered.  This  friction  allows  independent  lower- 
scale  control  of  horizontal  mass  divergence  and  momentum 
vorticity.* 

When  allowing  convectional  mixing,  the  heating,  particularly  the 
sensible  heating,  should  be  modified  accordingly. 

The  numerical  analogue  which  we  have  developed  is  also,  of  course, 

not  the  last  word.  We  have  tried  to  rigorously  couple  the  variables  by 

physical  relationship.  In  this  we  have  not  been  completely  successful, 

as  evidenced  by  some  double  differences  in  the  analogue.  The  selection 

of  levels  and  the  module  dimensional  constants  may  also  be  subjected  to 

further  statistical  refinements.  For  lengthy  evolutions,  the  leap-frog 

method  may  create  some  difficulties.  The  method  may  introduce  extraneous 

behavior  as  shown  by  linear  analysis.  A  non-linear  numerical  instability 

(8  ) 

may  also  arise;  however,  the  type  discovered  and  analyzed  by  Phillips 
in  1959  is  due  to  interactions  in  two  horizontal  dimensions.  The  use 
of  the  centered  single-difference  analogue,  Eq.  (83),  should  also  be 
considered . 

With  these  modifications,  it  should  be  possible  and  profitable  to 
carry  the  integration  through  a  complete  annual  cycle.  This  should 
be  possible  on  one  of  the  larger  available  computers,  in  the  order  of  a 
few  seconds  per  time  step.  The  evolution  can  again  be  started  with  a 
Standard  Atmosphere  at  rest,  or  with  balanced  zonal  motion.  The  heating 
should  be  introduced  in  steps  over  a  period  of  several  days,  in  order 
to  'void  exciting  a  48-hour  pole-to-pole  compression-gravity  oscillation. 

Such  an  experiment  would  yield  an  abundance  of  information  on  the 
primary  scale  (the  planetary  scale)  in  our  atmosphere.  It  would  indicate 
the  j.rculation  undisturbed  by  longitudinal  variations  and  dynamic- 
instability  breakdowns,  and  would  allow  a  separate  appreciation  of  these 
effects.  It  would  yield  more  understanding  of  the  cascade  of  energy 
downscale  from  the  primary  and  dictating  planetary  scale. 


* 


This  form  of  friction  with  experimental  findings  will  be  discussed 
in  a  forthcoming  report  by  the  author. 
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The  model  nay  also  be  used  to  investigate  a  number  of  geophysical 
influences  and  may  be  applied  to  arbitrary  planetary  atmospheres. 
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